---
title: 'Prevalence of trypanosomes and Sodalis in wild populations of tsetse flies:
  Impact of SIT programmes for tsetse eradication'
author: Mouhamadou M. Dieng, Mikhailou K. Dera, Percy Moyab, Gisele M. S. Ouedraogo,
  Guler Demirbas-Uzel, Fabian  Gstöttenmayer, Fernando C. Mulandane, Luis Neves,
  Sihle Mdluli, Jean - Baptiste Rayaisse, Soumaïla Pagabeleguem, Chantel J.
  de Beer, Andrew G. Parker,  Jan Van Den Abbeele, Robert L. Mach1, Marc J.
  B. Vreysen, and Adly M. M. Abd-Alla
date: "03/06/2021"
output: word_document
---

```{r}
setwd("C:/Users/abdallaa/OneDrive - IAEA/My_passport_6/Dieng_M/Sodalis_tryp_MS_20210603/Raw_data")

library(ggplot2)
library(lattice)
library(gcookbook)
library(ggfortify)
library(datasets)
library(MASS)
library(survival)
library(rmarkdown)
library(knitr)
library(coxme)
library(lme4)
library(nlme)
library(tidyverse)
library(gapminder)
library(rcompanion)
library(FSA)
library(stats)
library(RCA)
library(broom)
library(sp)
library(MuMIn)
library(ggpubr)
library(AICcmodavg)
library(car)
library(ggthemes)
```

## prepare the data 

```{r}
data=read.csv("rawdata_statistic_sod_tryp_bio2.csv")
str(data)
attach(data)
head(data)
data=na.omit(data)
data

data$Country=as.factor(data$Country)
data$Localisation=as.factor(data$Localisation)
data$Species=as.factor(data$Species)
data$Sex=as.factor(data$Sex)
str(data)
attach(data)
head(data)
data=na.omit(data)
data

```

## Statistics showen in the mansucript consequently

```{r}
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_Tspp) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_Tspp) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_Tspp) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_Tspp) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_Tspp) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model1)
Anova(model1)

Anova(model2)
#-----------------------------------------------------------------------
#impact of taxa and countries

# statistics in the manuscript
data_ga <- subset(data, Species=="Ga")
data_ga

model1<-glm((Prev_Tspp) ~ Country,data=data_ga, family=gaussian())
summary(model1)
Anova(model1)

#----------------------------------------------
data_gb <- subset(data, Species=="Gb")
data_gb

model1<-glm((Prev_Tspp) ~ Country,data=data_gb, family=gaussian())
summary(model1)
Anova(model1)

#-------------------------------------------
data_gff <- subset(data, Species=="Gff")
data_gff

model1<-glm((Prev_Tspp) ~ Country,data=data_gff, family=gaussian())
summary(model1)
Anova(model1)

#-------------------------------------------
data_gmm <- subset(data, Species=="Gmm")
data_gmm

model1<-glm((Prev_Tspp) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
Anova(model1)

#-----------------------------------------------
data_gp <- subset(data, Species=="Gp")
data_gp

model1<-glm((Prev_Tspp) ~ Country,data=data_gp, family=gaussian())
summary(model1)
Anova(model1)

#-----------------------------------------------
data_gpg <- subset(data, Species=="Gpg")
data_gpg

model1<-glm((Prev_Tspp) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
Anova(model1)

#-----------------------------------------------
data_gt <- subset(data, Species=="Gt")
data_gt

model1<-glm((Prev_Tspp) ~ Country,data=data_gt, family=gaussian())
summary(model1)
Anova(model1)

#------------------------------------------------------
# all species

model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)
Anova(model1)

```

## Statistics for Table 2

```{r}
#Glm Tspp per country
data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_Tspp) ~ Country,data=data, family=gaussian())
summary(model1)

#--------------------------------------------------------------
# for Sodalis

# model selection
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_Sod) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_Sod) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_Sod) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_Sod) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_Sod) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model3)
Anova(model3)

Anova(model1)

Anova(model2)

data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_Sod) ~ Country,data=data, family=gaussian())
summary(model1)

```

## Statistics for table 3

```{r}
# for trypanosome Tspp

data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_Tspp) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_Tspp) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_Tspp) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_Tspp) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_Tspp) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_Tspp) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_Tspp) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_Tspp) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_Tspp) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_Tspp) ~ Species,data=data, family=gaussian())
summary(model1)

#--------------------------------------------------------------
# for Sodalis

data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_Sod) ~ Species,data=data, family=gaussian())
summary(model1)
Anova(model1)

data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_Sod) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_Sod) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_Sod) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_Sod) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_Sod) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_Sod) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_Sod) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_Sod) ~ Species,data=data, family=gaussian())
summary(model1)

data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_Sod) ~ Species,data=data, family=gaussian())
summary(model1)

```

## Statistics for table 4

```{r}
# species with significant differences between countries

data_ga <- subset(data, Species=="Ga")
data_ga

model1<-glm((Prev_Tspp) ~ Country,data=data_ga, family=gaussian())
summary(model1)
Anova(model1)

model2<-glm((Prev_Sod) ~ Country,data=data_ga, family=gaussian())
summary(model2)
Anova(model2)
#-----------------------------------------------

data_gb <- subset(data, Species=="Gb")
data_gb

model1<-glm((Prev_Tspp) ~ Country,data=data_gb, family=gaussian())
summary(model1)
Anova(model1)

model2<-glm((Prev_Sod) ~ Country,data=data_gb, family=gaussian())
summary(model2)
Anova(model2)

#-------------------------------------------------------------------------
data_gff <- subset(data, Species=="Gff")
data_gff

model1<-glm((Prev_Tspp) ~ Country,data=data_gff, family=gaussian())
summary(model1)
Anova(model1)

model2<-glm((Prev_Sod) ~ Country,data=data_gff, family=gaussian())
summary(model2)
Anova(model2)

#-------------------------------------------------------------------------
data_gmm <- subset(data, Species=="Gmm")
data_gmm

model1<-glm((Prev_Tspp) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
Anova(model1)

model2<-glm((Prev_Sod) ~ Country,data=data_gmm, family=gaussian())
summary(model2)
Anova(model2)

data_gmm$Country <- relevel(data_gmm$Country, ref= "URT")
model2<-glm((Prev_Sod) ~ Country,data=data_gmm, family=gaussian())
summary(model2)
Anova(model2)

data_gmm$Country <- relevel(data_gmm$Country, ref= "ZAM")
model2<-glm((Prev_Sod) ~ Country,data=data_gmm, family=gaussian())
summary(model2)
Anova(model2)

data_gmm$Country <- relevel(data_gmm$Country, ref= "ZIM")
model2<-glm((Prev_Sod) ~ Country,data=data_gmm, family=gaussian())
summary(model2)
Anova(model2)
#----------------------------------------------------------------
#Gp

data_gp <- subset(data, Species=="Gp")
data_gp

model1<-glm((Prev_Tspp) ~ Country,data=data_gp, family=gaussian())
summary(model1)
Anova(model1)

model2<-glm((Prev_Sod) ~ Country,data=data_gp, family=gaussian())
summary(model2)
Anova(model2)

data_gp$Country <- relevel(data_gp$Country, ref= "KEN")
model2<-glm((Prev_Sod) ~ Country,data=data_gp, family=gaussian())
summary(model2)
Anova(model2)

data_gp$Country <- relevel(data_gp$Country, ref= "URT")
model2<-glm((Prev_Sod) ~ Country,data=data_gp, family=gaussian())
summary(model2)
Anova(model2)

data_gp$Country <- relevel(data_gp$Country, ref= "UGA")
model2<-glm((Prev_Sod) ~ Country,data=data_gp, family=gaussian())
summary(model2)
Anova(model2)

data_gp$Country <- relevel(data_gp$Country, ref= "ZAM")
model2<-glm((Prev_Sod) ~ Country,data=data_gp, family=gaussian())
summary(model2)
Anova(model2)

data_gp$Country <- relevel(data_gp$Country, ref= "ZIM")
model2<-glm((Prev_Sod) ~ Country,data=data_gp, family=gaussian())
summary(model2)
Anova(model2)

#------------------------------------------------------------------
#Gpg

data_gpg <- subset(data, Species=="Gpg")
data_gpg

model1<-glm((Prev_Tspp) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
Anova(model1)

model2<-glm((Prev_Sod) ~ Country,data=data_gpg, family=gaussian())
summary(model2)
Anova(model2)

data_gpg$Country <- relevel(data_gpg$Country, ref= "GUI")
model1<-glm((Prev_Tspp) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
Anova(model1)


data_gpg$Country <- relevel(data_gpg$Country, ref= "MLI")
model1<-glm((Prev_Tspp) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
Anova(model1)

data_gpg$Country <- relevel(data_gpg$Country, ref= "SEN")
model1<-glm((Prev_Tspp) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
Anova(model1)

#-------------------------------------------------------------
#Gt

data_gt <- subset(data, Species=="Gt")
data_gt

model1<-glm((Prev_Tspp) ~ Country,data=data_gt, family=gaussian())
summary(model1)
Anova(model1)

model2<-glm((Prev_Sod) ~ Country,data=data_gt, family=gaussian())
summary(model2)
Anova(model2)

```

## Selecting the GLM model for trypanosome species and mixed infection 

```{r}
#TC selection model
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_Tc) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_Tc) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_Tc) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_Tc) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_Tc) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model1)
Anova(model1)

#model2 is the best AICc=1084.04
Anova(model2)
#-------------------------------------------------------------------------
#Tv selection model
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_Tv) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_Tv) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_Tv) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_Tv) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_Tv) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model1)
Anova(model1)
#model1 is the best AICc= 1360.14     
Anova(model2)

#-------------------------------------------------------------------------
#Tz selection model
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_Tz) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_Tz) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_Tz) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_Tz) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_Tz) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model1)
Anova(model1)
#model1 is the best AICc= 1267.24
Anova(model2)

#-------------------------------------------------------------------------
#Tsg selection model
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_Tsg) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_Tsg) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_Tsg) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_Tsg) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_Tsg) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model1)
Anova(model1)
#model2 is the best AICc=1210.49
Anova(model2)

#-------------------------------------------------------------------------
#TcTv selection model
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_TcTv) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_TcTv) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_TcTv) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_TcTv) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_TcTv) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model1)
Anova(model1)
#model1 is the best AICc= 489.05
Anova(model2)


#-------------------------------------------------------------------------
#TcTz selection model
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_TcTz) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_TcTz) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_TcTz) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_TcTz) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_TcTz) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model1)
Anova(model1)
#model2 is the best AICc= 770.69       
Anova(model2)

#-------------------------------------------------------------------------
#TcTsg selection model
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_TcTsg) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_TcTsg) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_TcTsg) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_TcTsg) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_TcTsg) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model1)
Anova(model1)
#Model3 is the best AICc = 750.40       
Anova(model2)

#-------------------------------------------------------------------------
#TvTz selection model
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_TvTz) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_TvTz) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_TvTz) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_TvTz) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_TvTz) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model1)
Anova(model1)
#model1 is the best AICc = 1149.47       
Anova(model2)

#---------------------------------------------------------------------------
#TvTsg selection model
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_TvTsg) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_TvTsg) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_TvTsg) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_TvTsg) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_TvTsg) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model1)
Anova(model1)
#model4 is the best AICc = 1088.66
Anova(model2)

#----------------------------------------------------------------------------

#TzTsg selection model
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_TzTsg) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_TzTsg) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_TzTsg) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_TzTsg) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_TzTsg) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model1)
Anova(model1)
#model2 is the best AICc = 886.37       
Anova(model2)

#----------------------------------------------------------------------------

#TcTvTz selection model
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
model2<-glm((Prev_TcTvTz) ~ Species,data=data, family=gaussian())
model3<-glm((Prev_TcTvTz) ~ Country*Species,data=data, family=gaussian())
model4<-glm((Prev_TcTvTz) ~ Country+Species,data=data, family=gaussian())
model5<-glm((Prev_TcTvTz) ~ Country*Species+Localisation,data=data, family=gaussian())
model6<-glm((Prev_TcTvTz) ~ Country+Species+Localisation,data=data, family=gaussian())

#AICc(model1, model2, model3, model4)
model.set <- list(model1, model2, model3, model4, model5, model6)
model.names <- c("model1", "model2","model3", "model4", "model5", "model6")

aictab(model.set, modnames = model.names)
summary(model1)
Anova(model1)
#model2 is the best AICc = 265.72 
Anova(model2)

```

## Statistics for Supplementary table 3

```{r}
#======== Glm Tc per country
data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_Tc) ~ Country,data=data, family=gaussian())
summary(model1)

#======== Glm Tv per country
data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_Tv) ~ Country,data=data, family=gaussian())
summary(model1)

#======== Glm Tz per country
data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_Tz) ~ Country,data=data, family=gaussian())
summary(model1)

#======== Glm Tsg per country
data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_Tsg) ~ Country,data=data, family=gaussian())
summary(model1)

#======== Glm TcTv per country
data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_TcTv) ~ Country,data=data, family=gaussian())
summary(model1)

#======== Glm TcTz per country
data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)

#======== Glm TcTsg per country
data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_TcTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_TcTsg) ~ Country,data=data, family=gaussian())
summary(model1)

#======== Glm TvTz per country
data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_TvTz) ~ Country,data=data, family=gaussian())
summary(model1)

#======== Glm TvTsg per country
data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)

#======== Glm TzTsg per country
data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_TvTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_TzTsg) ~ Country,data=data, family=gaussian())
summary(model1)

#======== Glm TcTvTz per country
data$Country <- relevel(data$Country, ref= "BKF")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GHA")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "GUI")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ETH")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "KEN")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MLI")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "MOZ")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAI")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SEN")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SWA")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "URT")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "UGA")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "SAF")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZAM")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)
data$Country <- relevel(data$Country, ref= "ZIM")
model1<-glm((Prev_TcTvTz) ~ Country,data=data, family=gaussian())
summary(model1)

```

## Statistics for Supplementary table 4

```{r}
#====== Glm Tc per species
data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_Tc) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_Tc) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_Tc) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_Tc) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_Tc) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_Tc) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_Tc) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_Tc) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_Tc) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_Tc) ~ Species,data=data, family=gaussian())
summary(model1) 

#====== Glm Tv per species
data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_Tv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_Tv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_Tv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_Tv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_Tv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_Tv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_Tv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_Tv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_Tv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_Tv) ~ Species,data=data, family=gaussian())
summary(model1)

#====== glm Tz per species
data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_Tz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_Tz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_Tz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_Tz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_Tz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_Tz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_Tz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_Tz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_Tz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_Tz) ~ Species,data=data, family=gaussian())
summary(model1)

#======= Glm Tsg per species
data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_Tsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_Tsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_Tsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_Tsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_Tsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_Tsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_Tsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_Tsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_Tsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_Tsg) ~ Species,data=data, family=gaussian())
summary(model1)

#========= Glm TcTv per species
data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_TcTv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_TcTv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_TcTv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_TcTv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_Tsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_TcTv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_TcTv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_TcTv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_TcTv) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_TcTv) ~ Species,data=data, family=gaussian())
summary(model1)

#========= Glm TcTz per species
data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_TcTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_TcTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_TcTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_TcTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_TcTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_TcTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_TcTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_TcTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_TcTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_TcTz) ~ Species,data=data, family=gaussian())
summary(model1)

#======= Glm TcTsg per species
data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_TcTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_TcTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_TcTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_TcTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_TcTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_TcTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_TcTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_TcTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_TcTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_TcTsg) ~ Species,data=data, family=gaussian())
summary(model1)

#======== Glm TvTsg per species
data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_TvTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_TvTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_TvTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_TvTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_TvTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_TvTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_TvTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_TvTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_TvTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_TvTsg) ~ Species,data=data, family=gaussian())
summary(model1)

#========= Glm TvTz per species
data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_TvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_TvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_TvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_TvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_TvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_TvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_TvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_TvTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_TvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_TvTz) ~ Species,data=data, family=gaussian())
summary(model1)

#========= glm TzTsg per species
data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_TzTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_TzTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_TzTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_TzTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_TzTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_TzTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_TzTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_TzTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_TzTsg) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_TzTsg) ~ Species,data=data, family=gaussian())
summary(model1)

#======== Glm TcTvTz per species
data$Species <- relevel(data$Species, ref= "Ga")
model1<-glm((Prev_TcTvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gb")
model1<-glm((Prev_TcTvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gff")
model1<-glm((Prev_TcTvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmed")
model1<-glm((Prev_TcTvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmm")
model1<-glm((Prev_TcTvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gmsm")
model1<-glm((Prev_TcTvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gp")
model1<-glm((Prev_TcTvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpg")
model1<-glm((Prev_TcTvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gpp")
model1<-glm((Prev_TcTvTz) ~ Species,data=data, family=gaussian())
summary(model1)
data$Species <- relevel(data$Species, ref= "Gt")
model1<-glm((Prev_TcTvTz) ~ Species,data=data, family=gaussian())
summary(model1)

```

## Statistics for Supplementary table 5

```{r}
###For Ga
data_ga <- subset(data, Species=="Ga")
data_ga

##Tc
data_ga$Country <- relevel(data_ga$Country, ref= "SAF")
model1<-glm((Prev_Tc) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "MOZ")
model1<-glm((Prev_Tc) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "SWA")
model1<-glm((Prev_Tc) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "URT")
model1<-glm((Prev_Tc) ~ Country,data=data_ga, family=gaussian())
summary(model1)

##Tv
data_ga$Country <- relevel(data_ga$Country, ref= "SAF")
model1<-glm((Prev_Tv) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "MOZ")
model1<-glm((Prev_Tv) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "SWA")
model1<-glm((Prev_Tv) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "URT")
model1<-glm((Prev_Tv) ~ Country,data=data_ga, family=gaussian())
summary(model1)

##Tsg
data_ga$Country <- relevel(data_ga$Country, ref= "SAF")
model1<-glm((Prev_Tsg) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "MOZ")
model1<-glm((Prev_Tsg) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "SWA")
model1<-glm((Prev_Tsg) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "URT")
model1<-glm((Prev_Tsg) ~ Country,data=data_ga, family=gaussian())
summary(model1)

##TZ
data_ga$Country <- relevel(data_ga$Country, ref= "SAF")
model1<-glm((Prev_Tz) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "MOZ")
model1<-glm((Prev_Tz) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "SWA")
model1<-glm((Prev_Tz) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "URT")
model1<-glm((Prev_Tz) ~ Country,data=data_ga, family=gaussian())
summary(model1)

##TcTv
data_ga$Country <- relevel(data_ga$Country, ref= "SAF")
model1<-glm((Prev_TcTv) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "MOZ")
model1<-glm((Prev_TcTv) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "SWA")
model1<-glm((Prev_TcTv) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "URT")
model1<-glm((Prev_TcTv) ~ Country,data=data_ga, family=gaussian())
summary(model1)

##TcTsg
data_ga$Country <- relevel(data_ga$Country, ref= "SAF")
model1<-glm((Prev_TcTsg) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "MOZ")
model1<-glm((Prev_TcTsg) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "SWA")
model1<-glm((Prev_TcTsg) ~ Country,data=data_ga, family=gaussian())
summary(model1)
data_ga$Country <- relevel(data_ga$Country, ref= "URT")
model1<-glm((Prev_TcTsg) ~ Country,data=data_ga, family=gaussian())
summary(model1)

#----------------------------------------------------------------------------
###For Gb
data_gb <- subset(data, Species=="Gb")
data_gb


##Tc
data_gb$Country <- relevel(data_gb$Country, ref= "SAF")
model1<-glm((Prev_Tc) ~ Country,data=data_gb, family=gaussian())
summary(model1)

##Tv
data_gb$Country <- relevel(data_gb$Country, ref= "SAF")
model1<-glm((Prev_Tv) ~ Country,data=data_gb, family=gaussian())
summary(model1)

##Tsg
data_gb$Country <- relevel(data_gb$Country, ref= "SAF")
model1<-glm((Prev_Tsg) ~ Country,data=data_gb, family=gaussian())
summary(model1)

##Tz
data_gb$Country <- relevel(data_gb$Country, ref= "SAF")
model1<-glm((Prev_Tz) ~ Country,data=data_gb, family=gaussian())
summary(model1)

##TvTsg
data_gb$Country <- relevel(data_gb$Country, ref= "SAF")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gb, family=gaussian())
summary(model1)

#---------------------------------------------------------------------------
###For Gff
data_gff <- subset(data, Species=="Gff")
data_gff


##Tv
data_gff$Country <- relevel(data_gff$Country, ref= "KEN")
model1<-glm((Prev_Tv) ~ Country,data=data_gff, family=gaussian())
summary(model1)

##Tc
data_gff$Country <- relevel(data_gff$Country, ref= "KEN")
model1<-glm((Prev_Tc) ~ Country,data=data_gff, family=gaussian())
summary(model1)

##Tz
data_gff$Country <- relevel(data_gff$Country, ref= "KEN")
model1<-glm((Prev_Tz) ~ Country,data=data_gff, family=gaussian())
summary(model1)

##Tsg
data_gff$Country <- relevel(data_gff$Country, ref= "KEN")
model1<-glm((Prev_Tsg) ~ Country,data=data_gff, family=gaussian())
summary(model1)

##TcTv
data_gff$Country <- relevel(data_gff$Country, ref= "KEN")
model1<-glm((Prev_TcTv) ~ Country,data=data_gff, family=gaussian())
summary(model1)

##TcTz
data_gff$Country <- relevel(data_gff$Country, ref= "KEN")
model1<-glm((Prev_TcTz) ~ Country,data=data_gff, family=gaussian())
summary(model1)

##TcTsg
data_gff$Country <- relevel(data_gff$Country, ref= "KEN")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gff, family=gaussian())
summary(model1)

#----------------------------------------------------------------------
###For Gmm
data_gmm <- subset(data, Species=="Gmm")
data_gmm

##Tc
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZAM")
model1<-glm((Prev_Tc) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "KEN")
model1<-glm((Prev_Tc) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "URT")
model1<-glm((Prev_Tc) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZIM")
model1<-glm((Prev_Tc) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
##Tv
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZAM")
model1<-glm((Prev_Tv) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "KEN")
model1<-glm((Prev_Tv) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "URT")
model1<-glm((Prev_Tv) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZIM")
model1<-glm((Prev_Tv) ~ Country,data=data_gmm, family=gaussian())
summary(model1)

##Tsg
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZAM")
model1<-glm((Prev_Tsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "KEN")
model1<-glm((Prev_Tsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "URT")
model1<-glm((Prev_Tsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZIM")
model1<-glm((Prev_Tsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)

##Tz
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZAM")
model1<-glm((Prev_Tz) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "KEN")
model1<-glm((Prev_Tz) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "URT")
model1<-glm((Prev_Tz) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZIM")
model1<-glm((Prev_Tz) ~ Country,data=data_gmm, family=gaussian())
summary(model1)

##TcTv
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZAM")
model1<-glm((Prev_TcTv) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "KEN")
model1<-glm((Prev_TcTv) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "URT")
model1<-glm((Prev_TcTv) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZIM")
model1<-glm((Prev_TcTv) ~ Country,data=data_gmm, family=gaussian())
summary(model1)

##TcTsg
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZAM")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "KEN")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "URT")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZIM")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)

##TvTsg
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZAM")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "KEN")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "URT")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZIM")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)

##TzTsg
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZAM")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "KEN")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "URT")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)
data_gmm$Country <- relevel(data_gmm$Country, ref= "ZIM")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gmm, family=gaussian())
summary(model1)

#-------------------------------------------------------------------------------
###For Gp
data_gp <- subset(data, Species=="Gp")
data_gp

##Tc
data_gp$Country <- relevel(data_gp$Country, ref= "ETH")
model1<-glm((Prev_Tc) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "KEN")
model1<-glm((Prev_Tc) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "UGA")
model1<-glm((Prev_Tc) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "URT")
model1<-glm((Prev_Tc) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZIM")
model1<-glm((Prev_Tc) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZAM")
model1<-glm((Prev_Tc) ~ Country,data=data_gp, family=gaussian())
summary(model1)

##Tv
data_gp$Country <- relevel(data_gp$Country, ref= "ETH")
model1<-glm((Prev_Tv) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "KEN")
model1<-glm((Prev_Tv) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "UGA")
model1<-glm((Prev_Tv) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "URT")
model1<-glm((Prev_Tv) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZIM")
model1<-glm((Prev_Tv) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZAM")
model1<-glm((Prev_Tv) ~ Country,data=data_gp, family=gaussian())
summary(model1)

##Tsg
data_gp$Country <- relevel(data_gp$Country, ref= "ETH")
model1<-glm((Prev_Tsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "KEN")
model1<-glm((Prev_Tsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "UGA")
model1<-glm((Prev_Tsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "URT")
model1<-glm((Prev_Tsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZIM")
model1<-glm((Prev_Tsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZAM")
model1<-glm((Prev_Tsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)

##TZ
data_gp$Country <- relevel(data_gp$Country, ref= "ETH")
model1<-glm((Prev_Tz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "KEN")
model1<-glm((Prev_Tz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "UGA")
model1<-glm((Prev_Tz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "URT")
model1<-glm((Prev_Tz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZIM")
model1<-glm((Prev_Tz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZAM")
model1<-glm((Prev_Tz) ~ Country,data=data_gp, family=gaussian())
summary(model1)

##TcTv
data_gp$Country <- relevel(data_gp$Country, ref= "ETH")
model1<-glm((Prev_TcTv) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "KEN")
model1<-glm((Prev_TcTv) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "UGA")
model1<-glm((Prev_TcTv) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "URT")
model1<-glm((Prev_TcTv) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZIM")
model1<-glm((Prev_TcTv) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZAM")
model1<-glm((Prev_TcTv) ~ Country,data=data_gp, family=gaussian())
summary(model1)

##TcTZ
data_gp$Country <- relevel(data_gp$Country, ref= "ETH")
model1<-glm((Prev_TcTz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "KEN")
model1<-glm((Prev_TcTz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "UGA")
model1<-glm((Prev_TcTz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "URT")
model1<-glm((Prev_TcTz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZIM")
model1<-glm((Prev_TcTz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZAM")
model1<-glm((Prev_TcTz) ~ Country,data=data_gp, family=gaussian())
summary(model1)

##TvTZ
data_gp$Country <- relevel(data_gp$Country, ref= "ETH")
model1<-glm((Prev_TvTz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "KEN")
model1<-glm((Prev_TvTz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "UGA")
model1<-glm((Prev_TvTz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "URT")
model1<-glm((Prev_TvTz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZIM")
model1<-glm((Prev_TvTz) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZAM")
model1<-glm((Prev_TvTz) ~ Country,data=data_gp, family=gaussian())
summary(model1)

##TcTsg
data_gp$Country <- relevel(data_gp$Country, ref= "ETH")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "KEN")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "UGA")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "URT")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZIM")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZAM")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)

##TvTsg
data_gp$Country <- relevel(data_gp$Country, ref= "ETH")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "KEN")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "UGA")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "URT")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZIM")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZAM")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)

##TzTsg
data_gp$Country <- relevel(data_gp$Country, ref= "ETH")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "KEN")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "UGA")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "URT")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZIM")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)
data_gp$Country <- relevel(data_gp$Country, ref= "ZAM")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gp, family=gaussian())
summary(model1)

#---------------------------------------------------------------------
###For Gpg
data_gpg <- subset(data, Species=="Gpg")
data_gpg

##Tc
data_gpg$Country <- relevel(data_gpg$Country, ref= "BKF")
model1<-glm((Prev_Tc) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "GUI")
model1<-glm((Prev_Tc) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "MLI")
model1<-glm((Prev_Tc) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "SEN")
model1<-glm((Prev_Tc) ~ Country,data=data_gpg, family=gaussian())
summary(model1)

##Tv
data_gpg$Country <- relevel(data_gpg$Country, ref= "BKF")
model1<-glm((Prev_Tv) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "GUI")
model1<-glm((Prev_Tv) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "MLI")
model1<-glm((Prev_Tv) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "SEN")
model1<-glm((Prev_Tv) ~ Country,data=data_gpg, family=gaussian())
summary(model1)

##Tsg
data_gpg$Country <- relevel(data_gpg$Country, ref= "BKF")
model1<-glm((Prev_Tsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "GUI")
model1<-glm((Prev_Tsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "MLI")
model1<-glm((Prev_Tsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "SEN")
model1<-glm((Prev_Tsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)

##Tz
data_gpg$Country <- relevel(data_gpg$Country, ref= "BKF")
model1<-glm((Prev_Tz) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "GUI")
model1<-glm((Prev_Tz) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "MLI")
model1<-glm((Prev_Tz) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "SEN")
model1<-glm((Prev_Tz) ~ Country,data=data_gpg, family=gaussian())
summary(model1)

##TcTz
data_gpg$Country <- relevel(data_gpg$Country, ref= "BKF")
model1<-glm((Prev_TcTz) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "GUI")
model1<-glm((Prev_TcTz) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "MLI")
model1<-glm((Prev_TcTz) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "SEN")
model1<-glm((Prev_TcTz) ~ Country,data=data_gpg, family=gaussian())
summary(model1)

##TcTsg
data_gpg$Country <- relevel(data_gpg$Country, ref= "BKF")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "GUI")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "MLI")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "SEN")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)

##TvTz
data_gpg$Country <- relevel(data_gpg$Country, ref= "BKF")
model1<-glm((Prev_TvTz) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "GUI")
model1<-glm((Prev_TvTz) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "MLI")
model1<-glm((Prev_TvTz) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "SEN")
model1<-glm((Prev_TvTz) ~ Country,data=data_gpg, family=gaussian())
summary(model1)

##TvTsg
data_gpg$Country <- relevel(data_gpg$Country, ref= "BKF")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "GUI")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "MLI")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "SEN")
model1<-glm((Prev_TvTsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)

##TvTsg
data_gpg$Country <- relevel(data_gpg$Country, ref= "BKF")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "GUI")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "MLI")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)
data_gpg$Country <- relevel(data_gpg$Country, ref= "SEN")
model1<-glm((Prev_TzTsg) ~ Country,data=data_gpg, family=gaussian())
summary(model1)

#-------------------------------------------------------------------------------
###For Gt
data_gt <- subset(data, Species=="Gt")
data_gt

##Tc
data_gt$Country <- relevel(data_gt$Country, ref= "BKF")
model1<-glm((Prev_Tc) ~ Country,data=data_gt, family=gaussian())
summary(model1)
##TV
data_gt$Country <- relevel(data_gt$Country, ref= "BKF")
model1<-glm((Prev_Tv) ~ Country,data=data_gt, family=gaussian())
summary(model1)
##Tsg
data_gt$Country <- relevel(data_gt$Country, ref= "BKF")
model1<-glm((Prev_Tsg) ~ Country,data=data_gt, family=gaussian())
summary(model1)
##Tz
data_gt$Country <- relevel(data_gt$Country, ref= "BKF")
model1<-glm((Prev_Tz) ~ Country,data=data_gt, family=gaussian())
summary(model1)
##TcTV
data_gt$Country <- relevel(data_gt$Country, ref= "BKF")
model1<-glm((Prev_TcTv) ~ Country,data=data_gt, family=gaussian())
summary(model1)
##TcTz
data_gt$Country <- relevel(data_gt$Country, ref= "BKF")
model1<-glm((Prev_TcTz) ~ Country,data=data_gt, family=gaussian())
summary(model1)
##TcTsg
data_gt$Country <- relevel(data_gt$Country, ref= "BKF")
model1<-glm((Prev_TcTsg) ~ Country,data=data_gt, family=gaussian())
summary(model1)
##TvTz
data_gt$Country <- relevel(data_gt$Country, ref= "BKF")
model1<-glm((Prev_TvTz) ~ Country,data=data_gt, family=gaussian())
summary(model1)
##TcTvTz
data_gt$Country <- relevel(data_gt$Country, ref= "BKF")
model1<-glm((Prev_TcTvTz) ~ Country,data=data_gt, family=gaussian())
summary(model1)

```

## Preparation of Figure 4

```{r}
data_sodqpcr=read.csv("rawdata_fig4_dataqpcr.csv")
data_sodqpcr
str(data_sodqpcr)
attach(data_sodqpcr)
head(data_sodqpcr)
data_sodqpcr=na.omit(data_sodqpcr)

fig4.tiff<-ggplot(data_sodqpcr,aes(x=Sample ,y=log_copy, fill=Sample)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2)+ ylim(0, 8)
fig4.tiff

tiff("fig4.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(fig4.tiff+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = c(.95, .35),legend.justification = c("right", "top"))) + xlab(expression(bolditalic("Infection type"))) + ylab( expression (paste (bold("log10  "), bolditalic("Sodalis"), bold(" copy number"))))
dev.off()

model1<-glm(log_copy ~ Sample, data = data_sodqpcr)
summary(model1)
Anova(model1)
```

## Supplementary figure 2

```{r}
# supplementary figure 2A
Gmm <- subset(data_sodqpcr, Species=="G.m")
Gmm

sup_fig2a<-ggplot(Gmm,aes(x=Sample ,y=log_copy, fill=Sample)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2)+ ylim(0, 8)
sup_fig2a

tiff("sup_fig2a", width = 4, height = 4, units = 'in', res = 300)
plot(sup_fig2a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = c(.95, .35),legend.justification = c("right", "top"))) + xlab(expression(bolditalic("Infection type"))) + ylab( expression (paste (bold("log10  "), bolditalic("Sodalis"), bold(" copy number"))))
dev.off()

model1<-glm(log_copy ~ Sample, data = Gmm)
summary(model1)
Anova(model1)

#--------------------------------------------------------------------
# supplementary figure 2A
Gp <- subset(data_sodqpcr, Species=="G.p")
Gp

sup_fig2a<-ggplot(Gp,aes(x=Sample ,y=log_copy, fill=Sample)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2)+ ylim(0, 8)

model1<-glm(log_copy ~ Sample, data = Gp)
summary(model1)

tiff("sup_fig2a", width = 4, height = 4, units = 'in', res = 300)
plot(sup_fig2a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = c(.95, .35),legend.justification = c("right", "top"))) + xlab(expression(bolditalic("Infection type"))) + ylab( expression (paste (bold("log10  "), bolditalic("Sodalis"), bold(" copy number"))))
dev.off()

model1<-glm(log_copy ~ Sample, data = Gp)
summary(model1)
Anova(model1)

```

